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INTRODUCTION 


This report is the second of a pair of reports dealing 
with the formulative aspects of a capability for finite element 
thin shell instability analysis. The first report (Ref. 1) 
was devoted to development of the linear stiffness relation- 
ships for a doubly-curved triangular thin shell finite element. 
The present report describes the extension of these relation- 
ships to account for geometrically nonlinear phenomena, which 
encompass the instability effects of interest to this study, 
and the procedures used in calculation of a variety of non- 
linear solutions and instability load intensities. 

A large share of the literature on finite element analy- 
sis of the past decade has been devoted to the area of non- 
linear analysis. Due to this rapid and continuing rate of 
progress a number of state-of-the-art surveys have appeared. 
Martin (Ref. 2) examined this topic in 1969. Subsequently, 

Oden (Ref. 3) gave an assessment that covered a wider range 
of nonlinearities. Stricklin and co-workers have presented 
a number of reviews, one of the more recent being Ref. 4. 
Gallagher (Ref. 5) also has reviewed the aspects of solution 
procedures of geometrically nonlinear analysis. Because of 
this extensive available literature the present report does 
not attempt a review of the same material, although in the 
body of the report justification is generally given for the 
choice of particular procedures in preference to available 
alternatives . 

In addressing problems of geometrically nonlinear analy- 
sis a decision must be made relative to the coordinate system 
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to which all basic relationships are to be referenced. The 
Lagrangian, or fixed frame of reference is adopted in this 
report in preference to the Eulerian, or moving, coordinate 
system. The largest share of work in geometrically nonlinear 
finite element analysis has adopted the former, although it 
should be added that selection is most likely one of taste 
since insufficient evidence has been accumulated of the rela- 
tive efficiency of these alternatives in the total analysis 
scheme. Progress in the Eulerian approach is described by 
Yamada (Ref. 6). 

The present work is also restricted to small strain 
states, so as to preserve the applicability of the linear 
stress -strain law. 

The theoretical basis of the developments described here- 
in is the principle of stationary potential energy. This 
means that the fundamental relationships consist solely of 
the stress-strain law and the appropriate strain-displacement 
relationships. The transformation from continuum form to 
algebraic relationships is effected on the basis of element 
assumed displacement fields. This approach, of necessity, 
follows from the adoption of the same approach in the formula- 
tion of the linear stiffness relationships in Ref. 1. It is 
pertinent to note that alternative variational principles 
are not well developed in the area of geometrically nonlinear 
ahalysis, nor is it apparent that they would possess special 
advantages in finite element analyses of the subject class 
of problem. 
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The finite element formulation described in this report 
is not in strict adherence to the principle of stationary 
potential energy because of deficiencies in the satisfaction 
of conditions on interelement displacement continuity and on 
zero strain under rigid body motion. It is noted that their 
severity is reduced to acceptable proportions by means of 
interelement constraint conditions, and by the choice of 
the element displacement fields, and also that no new con- 
siderations in these respects are introduced by the extension 
to geometrically nonlinear analysis. 

The report is organized as follows. First, reviews are 
given of the types of geometrically nonlinear behavior likely 
to be sustained by thin shell structures and of the essential 
features of the linear stiffness formulation of the shell 
element treated herein. The latter enables one to read this 
report without recourse to Ref. 1. The formulation of both 
the element and global representations for the geometrically 
nonlinear finite element analysis of thin shells is then 
established without reference to a specific type of element, 
after which details are given of the formulation of the sub- 
ject triangular shell element. The nonlinear formulation is 
presented in both a rigorous, or "consistent" manner, in 
which the displacement field employed is that which was used 
in constructing the linear stiffness terms, and in an approxi- 
mate manner in which a simple linear displacement field is 
described over a quadrilateral region consisting of four tri- 
angles . 
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Two approaches to the solution of the global nonlinear 
algebraic equations, the Newton-Raphson and incremental 
methods, are discussed in some detail. The approach adopted, 
known as the modified incremental method, is one which com- 
bines one cycle of each of the above approaches within each 
load step. It is, of course, desirable to maintain control 
of the error produced in each step. For this purpose the 
same section of the report introduces a rational approach to 
error control, based on treatment of the analysis as an initial 
value problem. 

The two types of instability phenomena, limit points and 
bifurcat ional buckling, are studied in separate sections of 
the report. The limit point scheme is based on a superposition 
procedure applied within an increment of displacement, wherein 
all behavior is linearized. Bifurcation loads are determined 
through a Lagrangian interpolation scheme applied to the load 
versus tangent stiffness determinant solution points. A 
distinction is made between calculations for linear and non- 
linear states preceding bifurcation. 

The report concludes with a description of four numerical 
examples. The postbuckling response of a uniformly compressed 
flat plate with initial imperfections is studied for a variety 
of gridworks , solution algorithms, and support conditions. 
Analyses of a cylindrical shell segment under uniform radial 
load, a spherical cap under concentrated load at the crown and 
a clamped square plate under pressure are also described. Compari- 
sons are made in all cases with the results of available alterna- 
tive numerical or analytical solutions. 
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2. G EOMETRICALLY NONLINEAR BEHAVIOR OF SHELLS 

An understanding of the types of geometrically nonlinear 
response likely to be sustained by thin shell structures is 
essential to the construction of basic formulations and 
algorithms which allow such response to be predicted analyt- 
ically through the finite element method. 

The behavior of all structures at low load levels can 
be classified in one of three ways (Fig. 1). 

i) Linear load-deflection relationship 

ii) Nonlinear weakening 

iii) Nonlinear stiffening 

Common examples of weakening structures are shallow trusses, 
arches, or shells subjected to inwardly directed loads. The 
same structures are stiffening if the loads are directed out- 
wards. Examples of both types of structures are shown in 
Fig . 1 . 

The work described in this report is mainly concerned 
with nonlinear weakening situations and to describe in more 
detail the behavior to be anticipated, Fig. 2 has been drawn 
The solid line (OB) applies to "perfect" structures and, for 
instability problems, represents the case in which the struc- 
ture first displaces along the OA portion of the path (the 
fundamental path ) and bifurcates (or branches) at the point A 
Thus, uniqueness in the definition of the path is lacking at 
the point A. It should also be noted, in passing, that the 
assumption is often made in analytical approaches that the 
fundamental path is linear in the range OA. 
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The post -buckling path may rise (path AC, Fig. 2a) , or 
it may descend (path AD, Fig. 2b). The slope of the post- 
buckling path of the perfect structure gives an indication 
of the behavior to be expected of real structures, which in- 
evitably possess fabricat ional imperfections. For structures 
with fabricational imperfections the load-displacement be- 
havior in the presence of destabilizing phenomena follows 
the paths indicated by dotted lines. In such cases there is 
no bifurcation phenomenon. Structures with a rising post- 
buckling path as in Fig. 2a will have strength exceeding the 
bifurcation load. Imperfect structures with a descending 
post-buckling path in the perfect state (Fig. 2b) will not 
achieve strengths as high as the bifurcation load unless the 
load-displacement path again rises at larger displacements. 
Such structures, under the appropriate load condition, are 
therefore imperfection sensitive and the maximum load at- 
tained (Point E) is termed the limit point . The instability 
phenomenon which occurs at this point is known as snap- 
through, in identification with the physical circumstance 
of structures in this mode of load-displacement response. 

A non-bifurcating load-displacement behavior of the limit 
point (or snap-through) type may also occur for a structure 
devoid of imperfections. It should be cautioned that the 
realm of possibilities of structural instability phenomena 
iSs not exhausted by the circumstances defined in Fig. 2. 
Nevertheless, many common structures behave in this manner; 
some familiar examples are portrayed in Fig. 3. 
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3. REVIEW OF BASIC ELEMENT REPRESENTATION 

A brief review is given in this section of the features 
of the triangular thin shell element formulation for linear 
analysis. The review serves to define basic data and 
theoretical considerations needed for the extension of this 
formulation to accommodate geometrically nonlinear behavior. 
The three principal aspects of the linear formulation are 
Ca) the mode of geometric description, (b) the adopted strain- 
displacement equations, and (c) the selection of displacement 
functions. Also special considerations arise in the treatment 
of displacement incompatibilities along the element juncture 
lines. Each of these items is discussed in the following. 

The basic geometry of the triangular shell element is 
shown in Fig. 4. Points are located on the shell surface in 
the orthogonal curvilinear system a-3. The user provides 
nodal values for the thtee curvatures as well as the shell 
thickness t. These variables are interpolated linearly over 
the area of the element. 

The strain-displacement equations given by Koiter (Ref. 

7) are used in constructing the element potential energy for 
linear analysts. Koiter' s relationships are described as 
"consistent", i.e., terms of the order of the ratio of the 
membrane strain to the radius of curvature are retained or 
discarded in a consistent manner. Also, these relationships 
evidence no strain under rigid body motion. 

In discretizing the potential energy expression formed 
with use of the strain-displacement equations, the same dis- 
placement functions are chosen for u, v and w. Each is in 
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the form of a complete cubic polynomial, of 10 terms, for 
which the corresponding displacement fields are 


u = |_Nj{u} , v = j_N_j {v}, w = lNj{w} 


CD 


where l Nj = lNj .... N 4 J 


fu > ’ L u l CffJ (It) 


36' j u 2 


r3Us T 

^ 38^3 U ' 


l 4J 


( 2 ) 


and {v} and {w} are similarly defined. The subscript 4 refers 
to the centroidal node point. Details of these functions and 
the reasons for their choice are discussed in Ref. 1. As 
previously noted, neither of the two factors prominent in 
the choice of displacement fields for thin shell elements, 
the condition of zero strain energy under rigid body motion 
of the element nor that of continuity of angular displacement 
across element boundaries, are satisfied by the chosen func- 
tions, Eq. (1). The use of cubic functions for in-surface 
displacement components (u and v) represents a move in the 
direction of the implicit satisfaction of the former condi- 
tion. 

The failure to meet conditions on the interelement con- 
tinuity of angular displacement has more serious consequences 
on solution accuracy and for this reason special steps are 
taken to "restore" this continuity. The approach taken, 
and indeed the use of a cubic polynomial in representation 
of the flat triangular element, was introduced by Anderheggen 
(Ref. 34) and Harvey and Kelsey (Ref. 8). If A and B are 
neighboring elements and the subscript n denotes the normal 
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direction, the relative angular displacement of adjacent 
edges at their midpoint i can be set to zero to ensure inter- 
element continuity. 


A B 

, 9W-. _ ,3w. ,3Wv 


0 


( 3 ) 


Now, the differentiation of the displacement field w yields 
the angular displacement in terms of the joint displacements, 
so that a constraint equation of the following form can be 
written 

{c a-b >T{ a} = 0 (4) 


The above process is applied to each of the m element 
boundaries and for a system with n degrees-of-freedom the 
following m by n set of equations is obtained 

[C] {A} = 0 (5) 

Eq. (5) is supplemental to the basic force -displacement 
relationships of the stiffness analysis. The manner in which 
such constraint equations are used to modify the basic force- 
displacement relationships is taken up at the close of the 
next section. 

4 • FORMULATIO N OF N ONLINEAR P ROBLEM 

The derivation given in this section is based on the 
principle of stationary potential energy. The potential 
energy (lip) of a shell element can be written in the follow- 
ing general form 

IL = 1 / L ej [D] { e) dA + V 
p * A 


( 6 ) 
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where A is the surface area of the shell 

l_ej is the row vector of relevant strain components 
({e} = |_ej T ) 

[D] is the matrix of constants in the relationships 

between the strains and the corresponding stresses, 
in the form {a} = [D] {e}. 

V is the potential of the applied loads. 

For simplicity, a more explicit statement of the potential 
of the applied loads is deferred until Eq. (6) has been dis- 
cretized and developed in greater detail. 

The present work is concerned with geometric nonlinear- 
ities, whose representation is embodied in the strain- 
displacement equations. Such equations can be expressed in 
the following general form 

(e) = {e L } + { e N } (7) 


where {e^} gives the linear portion of the strain-displacement 
equations in terms of displacements and {e N } gives the non- 
linear portion. The present work adopts a Lagrangian frame 
of reference for the construction of these relationships, 
i.e., all quantities are referenced to the axes of the un- 
deformed state. Thus, for the illustrative example of an 
axial member (Fig. 5) 

L _ du 
e “ 


and 


e N _ l du. 2 1 r dv. 2 1 r dw. 

e ' 1 c ax } I C 3P W 


The present work also assumes that the direct and shear 
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strains are small, so that the first term in the expression 
N 

for e is discarded, as are similar terms in the expressions 
for {e } for shell behavior. 

By substitution of Eq. (7) into Eq. (6) there is 
obtained 

n p = \ J A [ L £ j IdH £L } ♦ l £ J lDKe N > 

+ L_e^j[D]{e L > + l £ ^ [D]{e N }]dA + V (7a) 

This expression will next be discretized with use of the 
same displacement fields chosen for the linear formulation, 

Eqs . (1). For conciseness the three components u, v and w 
are grouped under the symbol A^ and Eq. (1) can be restated 
as 

A = [N] { A } (8) 

where {A} is a 30x1 vector which lists the joint displace- 
ments in the order defined below Eq. Cl)* The 3x30 matrix 
[N] is similarly arrayed. By performance of the operations 
dictated by the strain-displacement relationships the dis- 
cretized form of the strains is given by 


{e L } = 

[B L ]{A> 

(9a) 

{e N > = 

[b n ]{a> 

(9b) 


The terms of [B^] are constants and functions of the surface 

N 

coordinates and the terms of [B ] are functions of the dis- 
placements {A) as well as of constants and the surface co- 
ordinates . 
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After substitution of Eqs. (9a) and (9b) into Eq. 

(7a) and performance of the indicated integration the finite 

element (discretized) potential energy II is obtained. In 

~P 

indicial notation 


II ~ k A. A. + -jr k.. A. A- A. + * k 9 A. A. A. A 0 

,p 2 -o ij -1 - 3 6 l ijk -I -J -k IT 2 ijki -1 -} -k 


+ V 


(10a) 


in which k , k, , and k 9 are second, third and fourth 
°ij 1 ijk z ijk£ 

order tensors and are fixed values for a given element, and 
V is the potential of any applied loads on the surface of 
the element. The subscript tilde (~) denotes that these are 
element quantities. Alternatively, in matrix notation 

5p ■ [k o ]{A> * [k 1 (A)]{iJ ♦ r | L - J Ck 2 CA 2 ) I {A> ♦ V 

(10b) 

2 

The matrices [k Q ] , [k^(A)] and [k 2 (A )] are the linear and 

first -order and second order geometric element stiffness 

? 

matrices, respectively. The symbols (A) and (A ) are inserted 
adjacent to k^ and k 2 to emphasize that terms of these matrices 
are, respectively, linear and quadratic functions of the dis- 
placements {A}. Elsewhere in this report, with certain excep- 
tions, these matrices will simply be denoted as [k^] and [k 2 ] . 
It is also noted that other reports of this research (Refs. 9 
and 10) have designated these matrices as [n^] and [n 2 ] . To 
avoid confusion with the designation of shape functions given 
in this and the companion report, however, t.hese matrices are 
designated as [k^] and [k 2 ] herein. 
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For simplicity consideration is given only to concentrated 
applied loads acting at the degrees -of-freedom A^. Also, 
for the conditions of the present work it is assumed that a 
fixed distribution of normalized loads P\ prevails at all in- 
tensities of load n so that the loads are alternatively writ- 
ten as P. a n The potential of the applied loads is there- 

fore V = X A^, where the summation convention of indicial 
notation prevails. The potential energy of the complete struc- 
ture is evaluated by summing the contributions of all elements 
and of V. In indicial notation 


n = =- K A. A. 

p 2 ° i j l 2 


+ 4- K-. A. A. A. 
5 1 ijk 1 J k 


IT K 2 ijk /i A j A k i li ' n Vi 


(11a) 


Alternatively, in matrix notation 

n P = r LAJ [* 0 Ha> + £- lAj [k 1 (a)]{a> + t ^ lAJ [k 2 (a 2 )]{a> 


- n LAj {F} ( 11 b) 

The matrices [ K Q ] , [K^(A)] and [K 2 (A 2 )] are the global linear 
and first and second order geometric stiffness matrices and 
again the dependence of the latter two on {A} has been 
emphasized. 

By taking the first variation of 11^ with respect to each 
degree-of-freedom (A^) in turn, and setting each to zero in 
accordance with the stationary property of potential energy, 
and noting that the first three terms on the right side are 
quadratic, cubic, and quart ic functions, respectively, of the 
degrees of freedom, one obtains the nonlinear algebraic equa- 
tions of equilibrium. In indicial notation 
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K o A i * \ K i Mk * T K ? Mk A i - 1 F i ■ 0 (12a) 

°ij 3 i I ijk 3 K a ^ijkl 3 K 1 • 1 


or, in matrix notation 

[K q ] { A } + j [K 1 3 {A} + j [K 2 ]{A> - n {P} = 0 (12b) 


It is to be observed that special care must be taken in the 

construction of the [K^j and [K 2 3 matrices so that the simple 

d L^J 

mode of differentiation indicated above (e.g. aTA} ( ~6“ [K l ]{A}) 
“ \ l K ].]{A}) holds true. This point is elaborated upon in 
the next section. 

No explicit consideration has been given in the foregoing 
to the constraint conditions (Eq . (5)) which are necessary for 
the restoration of interelement displacement continuity. The 
next section will show that these constraints do not influence 
the nonlinear matrices [K^ ] and [K 2 ] . Thus, only [K q ] is in- 
fluenced and it is assumed that the full linear portion of the 
equilibrium equations are in fact of the form 


L c °J 

5. FORMULATION OF ELEMENT NONLINEAR MATRICES 


a. Consistent -Triangle 

This section is devoted to a detailed description of the 
formulation of the nonlinear matrices [K^] and [K 2 ] , which 
were introduced in the previous section, for the specific 
case of the subject triangular thin shell element. A cor- 
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respondingly detailed description of the formulation of the 
linear stiffness matrix [k Q ] is given in Ref. 1. The present 
formulation adheres to the notion of M consistency M , wherein 
the displacement functions employed in the construction of 
[k^] and j^] are identically those employed for [k Q ] . This 
results in a high cost of formulation for [k^] and [k 2 ] , so 
in the latter portion of this section a more economical, in- 
consistent, formulation is described. 

To establish the desired formulations, it is appropriate 
to compare terms in Eqs. (7a) and (10b). On this basis it is 
found that (the tilde underline is neglected, since this 
clearly refers to the individual element) 

5- LAJ [k 1 ]{4} - / Ci. E ^ [D]{e L ) * L E^J [D]{e N })dA (13) 
A 

and 

rl 1 -*- 1 (k 2 ](A) - ; (_e N j [D](e N }dA (14) 

A 

In order to proceed further the specific form of the 
relevant strain-displacement equations is needed. As noted 
in Section 2 the adopted linear components of these relation- 
ships are those due to Koiter (Ref. 7) . The corresponding 
nonlinear terms have been derived by Mushtari (Ref. 11). The 
strain components of thin shell behavior are the direct 
strains of membrane action e ^ and <f> and for flexure they 

are the curvatures k^ , and t. Mushtari has pointed out 
that the neglect of nonlinear terms in the curvature expres- 
sions is justified in problems of "medium" bending, i.e.. 
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where the normal displacements are of the same order of magni- 
tude as the thickness of the shell but much smaller than the 
other dimensions. Adopting this assumption, only the membrane 
strain-displacement equations are influenced. Hence 


1 3ij v dA w 1_ r 3w> 

e l ~ A 3a AF 30 " 2A 2 l 3a J 


_ 1 3v u 3B w 1_ f 3ws 

e 2 B 33 AB 3a " FJ 2fi 2 l 36 J 


(15) 


, _ 1 3v 1 3u u 3B v 3B 2w 1 r 3w. r 3w«. 

* A 3a F 36 ~ AF 3a " AF 3a T IAF ( -3a J L 36 J 


The nonlinear terms are underlined for clarity. A and B are 
the metric coefficients of the shell surface and R^ , R 2 and T 
are the radii of curvature of the undeformed shell in the 
a and 6 directions. 

In view of Eqs. (15), we have 

L 13u x v 3A w 

' R i 

___ _ w 

2 B 38 T AB 3a ' 


A 3a AB 36 


e N s _L riwq 

1 2A 2 ^ 


1 3v . u 3B 


N 

e 2 


2B 


o 


1 3v . 1 3u 
A 3a B 36 


u 3B _ v 3B _ 2w 
AF 3a AB 3a ~T 


( 16 ) 


.N 1 r 3w, ,3w. 
* “ 2AB C 3a H 36' ) 
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By substitution of the shape functions (Eq. (1)) into 
the above, expressions of the following form are obtained 
(consistent with the symbolism of Eqs. (9a) and (9b)) 


:7 » | B~ B i 

1 L uv W J 


U vl 



I 

L w . 

1 



N 1 . 
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(18) 


L N L N 

and similar expressions for e 2 » ^ an<1 $ • 

When expressions (17) and (18) are substituted in Eq. (13) 

terms of the following form are encountered. 



(19) 


The summation of these terms can be designated in indicial 

form as k, A.A-A V . The calculation of k, , once and for 
1 ij k 1 3 K 1 i j k 

all, would seem to be an attractive approach. Unfortunately, 

excessive computer storage required for the tensor prohibits 

this. Because of this the structure of Eq. (19) is reduced 

to quadratic form so that the kernel [k^] is a function of 

the displacements. This is done in the following manner. 


L N L _ 1 r L, N, T N a L , N. T N . N , N. T L N r N. T L 
E i E j E k " F l e i< e j> E k + E i < e k> E j + E j (E k } E i + e j^ e k> e i 


N f L, T N a N , N. T L, 
E k (E i } E j + E k (E j } E i^ 


( 20 ) 
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Typically, by invoking the nomenclature of Eqs. (17) and 
(18), the first term on the right side of Eq. (20) can be 
written as 
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which is the desired quadratic form. This procedure of trans- 
forming the product of three strain quantities into quadratic 
form with a nonlinear kernel is applied in the same manner 
to develop all terms of [kjJ. A listing of all of the result- 
ing terms can be found in Ref. 12. Clearly, the integration 
of these terms cannot be performed explicitly, so that a 
numerical integration, must be adopted. A 4 x 4 Gaussian 
scheme is employed, as in the case of [k Q ] in Ref. 1. 

The number of matrix operations required can be reduced 
considerably in two ways. First, by taking advantage of sym- 
metry, and secondly by noting that many of the matrices 
present can be obtained by transposing others. 

The same approach as that delineated above is used in 
calculating [k 2 ] . In this case products of four strains are 
converted to a quadratic form. The basic components of the 
[k 2 J matrix have the following structure 
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In indicial notation this would be k- A i A i A k A £* The equa- 

ijkA J 

tion corresponding to Eq. (20) for four strains is 
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The first term on the right side is rewritten as 


(23) 



which again is the desired form. By substituting terms of 
this sort for each product of two nonlinear strains in the 
energy relationship (Eq. (14)), the [k 2 ] matrix is obtained. 
Details of the resulting matrix are presented in Ref. 12. 

It is recalled (see Ref. 1 and Sect. 4 of this report) 
that the linear stiffness matrix that is relevant to the 
present study contains three rows and columns that correspond 
to constraint equations to enforce interelement continuity. 


L 
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The coefficients in these equations are constants that depend 
only on the element shape. Since the inplane displacements 
are considered to be negligible, the change in element shape 
is discounted. Consequently the constraint equations need no 
modification during a nonlinear analysis. 

b. Inconsistent - Quadrilateral 

Considerable savings in data preparation and computer 
execution times can be realized by the use of quadrilateral 
elements as compared to triangular elements. In the present 
work four triangles are combined to form a quadrilateral 
(Fig. 6). Constant curvature is assumed for the new element. 
Internal degrees of freedom are eliminated from the linear 
stiffness matrix by the use of static condensation. For non- 
linear analysis the geometric stiffness matrices are computed 
directly for a quadrilateral without cons iderat ion of the 
component triangles as delineated in this section. 

The finite element method, when based in variational 
principles, requires interelement continuity of derivatives 
up to one order lower than appears in the associated function- 
al, or energy. For the case of plate bending or shell struc- 
tures continuity of slopes is required of the displacement 
fields. This conclusion is valid for linear terms. With 
reference to the strain displacement relationships the high- 
est derivative to appear in the nonlinear terms is the first. 
Thus, only continuity of the displacements themselves is 
demanded by the variational principle. This opens up the 
possibility of using a simplified field for the nonlinear 
terms as compared to the complex field used for linear terms. 
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This is the approach adopted herein for the quadrilateral 
element . 

The approach introduced by Mallet and Marcal (Ref. 14) 
is modified to accommodate the curvature of an element. The 
strain displacement relationships of Eq. (15) are simplified 
to read 
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These can be written concisely as (Ref. 13) (with <p = e^) 

e i " L L l_l td} + 7 L d j [H i ] {d> (15b) 
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By introducing Eq. (15a) into the energy expressions, 

Eqs. (13) and (14), "core" forms of the first-and second- 
order geometric stiffness matrices can be established. 

A A 

These are denoted by [k^j and [k 2 ] and are given by 

[k x ] = / t C ± . ({L.} |_dj [Hj ] + |_dj {L i )[H j ] 

A 

+ [Hj] {d}j^LjJ) dA (25a) 

fk 2 ] - / t C.J C[H.]{d} L d Jf H j] * ?L d J IHjHdJtHiJJdA (25b) 

A 

where C^- is the material compliance and t is the thickness 
of the shell. These matrices, which are independent of the 
displacement field postulated for the element, are detailed 
in Ref. 12. 

In the present development a bilinear assumption is made 
for each of the three components of displacement, i.e. 

0 = 8^^ + a 2 a + a 3 & + a 4 a3 

v = a 5 + a 6 a + a 7 $ + a g cx3 (26) 

w ^ pet ^ a^^ 3 ^ a^af? 

By appropriate differentiation, vector {d} can be expressed 
in terms of the coefficients {a}. Thus 

{d} = [G] {a} (27) 

It is also necessary to express the generalized coefficients 
{a} in terms of the nodal displacements {A}. This is done 
by evaluating the polynomials of Eq. (26) at the node points. 
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{a} = [ J] {A} (28) 

The ingredients for forming the geometric stiffness 
matrices with respect to the nodal displacements are now 
ready. Thus 

[k,] = [J] T J [G] T [L][G]dA [J] 

1 A 1 

- [J] T [kjHJ] (29a) 

and 

[k 2 ] = [J] T [k 2 ][J] (29b) 

The necessary integration operations for [k^] and [k 2 ] can be 
performed explicitly. This has been done in the present case, 
with results given in detail in Ref. 12. Because of this 
ability for calculation of explicit terms the calculation of 
the inconsistent matrices for the quadrilateral is far more 
rapid than the consistent formulation calculation for the 
triangle . 

6. SOLUTION OF NONLINEAR EQUATIONS 
a. Newton-Raphson Method 

As noted in the Introduction, algorithms for the solution 
of the nonlinear algebraic equations representing geometrical- 
ly nonlinear behavior (e.g., Eq. (12b)) have attracted a 
great deal of attention from researchers. No one algorithm 
has gained universal preference over the others, mainly be- 
cause each method is effective for particular situations but 
inefficient or invalid for different situations. The approach 
adopted in the present work is therefore a combination of two 
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basic schemes, the Newton -Raphson and incremental methods 
respectively, and each of these is outlined in the present 
section. Another view of pertinent algorithms is that the 
analysis can be treated as an initial value problem. This 
view gives a different perspective on the methods adopted 
and also is a basis of a procedure for rational selection 
of load increment selection. Hence, the initial value analy- 
sis concept is also discussed herein. A concluding portion 
of the section ties together the overall strategy for solu- 
tion of the pertinent nonlinear equations. 

To describe the Newton-Raphson approach it is convenient 
to rearrange Eq. (12b) as follows 

{£} - [k Q ] { A} + \ [kj] {A} + | [k 2 ] { A> - n (F) - o (30) 

where (f) is an imbalance of load. A nonlinear analysis must 
proceed on the basis of an estimate of {A} and improve on the 
estimate in succeeding iterations. By expanding in a Taylor 
series about an approximate {f^}, after i iterations 


8{f.} 


{f..,} = { f - } + [ — 1 — ] {dA} + higher-order terms in {dA}(31) 
1 + JL 2. d A 


where {dA} = {A^ +1 } " C } is the change of displacements. 

With neglect of higher-order terms and with {f^ +1 > set equal 
to zero, we have 

<W ■ <*i> • ff i> ( 32a > 


or using the nomenclature of equation 30 


-1 


{1 W ■ { V ■ [ WK 2 1 {f i } 


(32b) 
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The iteration continues until {A^ + j> is arbitrarily close 
to {A i >. 

It is apparent that the "inversion" of a large matrix 
is called for at each iteration. To circumvent this problem 
the matrix of stiffness coefficients can be held fixed for 
several iterations so that modifications only to the residual 
load vector {f} are required. This approach is usually 
called the modified Newton -Raphson method. Its relationship 
to the standard Newton-Raphson method for a one-degree -of - 
freedom problem is shown graphically in Figure 7. 

It is obvious from Figure 7 that the convergence charac- 
teristics of the modified approach are inferior to the standard 
method but, because each iteration is cheaper in computational 
effort, the overall cost could well be less. Two factors 
have a bearing on this outcome. First is the number of equa- 
tions in the system. For a large order system the expense of 
repeated 'inversion* becomes very high. On the other hand, 
because of the increased number of iterations for the modified 
approach, the geometric stiffness matrices have to be evaluated 
a greater number of times. For simple (usually inconsistent) 
formulations these element form times are insignificant, but 
for complex derivations the times can compound to be greater 
than the time needed for solving the linear equations at each 
iteration. It would therefore seem that the standard Newton- 
Raphson method should be used for an idealization using com- 
plex elements where there is a corresponding reduction in the 
number of freedoms. The modified approach would be more com- 
petitive for large problems using simple elements. 
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The main disadvantage of the Newton -Raphson method is 
that it cannot be used in problems for which the solution is 
path dependent such as plasticity. In that particular in- 
stance, elastic unloading could occur during the iterative 
process which was not experienced by the actual structure. 
Thus, this method is not suitable for a general purpose 
program that accommodates geometric and material nonlineari- 
ties. Another shortcoming of the method in the form described 
above is that it does not converge at points of instability. 
This is so because the tangential stiffness matrix is singular 
at such points. Thurston (Ref. 15) resolves this difficulty 
by use of higher order terms in the Taylor series expansion 
upon which the method is based (Eq. 31). A more popular ap- 
proach in finite element analysis has been to prescribe dis- 
placement increments, rather than load increments, at or near 
critical points. A form of this approach has been adopted in 
the present \vork and is described in Section 7. 

Although in theory it is possible to solve the nonlinear 
system at any load level without previous information, care 
should be taken in selecting a starting point for the itera- 
tive process. Divergence is possible if the starting point 
is not sufficiently close to the solution. To overcome this 
problem it is recommended that solutions be achieved at lower 
load levels so that a reliable starting value can be safely 
extrapolated, 
b. Incremental Method 


In this method a solution {A + dA} is sought for Eq. (12b) 
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close to a known solution {A}. Using a Taylor series expan- 
sion for [K^ ] and [K 2 3 we get 

[1^ (A+dA)] {A+dA} « [ [K^ (A) ] + ^ [K 1 (A)]{dA} + ...]{A+dA} 

= [ (K^ (A) ] {A} + 2 [K 1 (A)]{dA> + (33) 

and 

[K 2 (A+dA) 2 ] { A+dA} » [K 2 (A 2 )]{A} + 3 [K 2 (A) 2 ]{dA> + ... (34) 

Neglecting terms of higher order than those shown explicitly 
on the right hand side, substituting these expressions into 
Eq. (12b), and collecting terms 

[[K 0 ] + [K i 3 + [K 2 ]](dA} = (n+dn) { p> - [k q 3 { a} 

- \ [K 1 ]{A> - | [K 2 ]{A> = dr\{¥} + {f} (35) 

By assuming that equilibrium is satisfied exactly at a load 
level denoted by n this equation can be simplified by sub- 
stituting { f > = 0 

C [K 0 3 + [K x 3 + [K 2 3]{dA> = dn {F} 

or 

[K t 1 (dA) = dn (F) (36) 

This is the basic equation of the incremental method. 
Thus at each load level one ’inversion' of the stiffness 
matrix is required, but the displacements are directly avail- 
able without any iteration. This approach is ideally suited 
for tracing the load-deflection behavior of a structure, 
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especially in the presence of material nonlinearities where 
small load increments are required to follow the spread of 
plastic zones. 

The main inaccuracy inherent in the incremental method 
is the gradual divergence from the true solution caused by the 
piecewise linearizations of a curve. This is shown graphical- 
ly in Figure 8. This error can be reduced, but not eliminated 
completely, by using smaller load steps. Another effective 
device for keeping the divergence to a minimum is to add the 
out of balance force vector ^f][ of the previous load level 
to the right hand side at the present load level. When this 
operation is performed on Eq. (36) one arrives back at Eq. 
(35). The algorithm indicated by Eq. (35) is called the 
modified incremental method. 

Under close scrutiny it is apparent that Eq. (35) and 
Eq. (32b) of the Newton-Raphson method are almost identical. 
The only difference being that in the incremental method 
{dA} is solved directly without iteration. For this reason 
the modified incremental approach can be thought of as one 
cycle of the Newton-Raphson procedure. 

c. Initial Value Approach for Error Estimate 

The geometrically nonlinear problem can be formulated 
as an initial value problem. When this is done it is pos- 
sible to identify a variety of established solution tech- 
niques which can be applied to the conditions under study. 
Although these alternatives are not adopted herein it is 
useful to examine the basic statement of the initial value 
format. This format sheds some light on the methods outlined 
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above. Also, it furnishes the basis for a rational scheme 
of increment step selection. 

Eq. (36) can be written in the limit, as dn goes to 
zero, as 

[[K 0 J + [Kj] + EK 2 )]{A}={F> (37) 

where the dot indicates differentiation with respect to the 
load level n. The simplest method of solving this equation 
is Euler's one step method. This yields Eq. (36) of the in- 
cremental method, as is to be expected since Euler's formula 
and Eq. (36) are both derived by expanding about a known solu- 
tion point by a Taylor series expansion, quadratic and higher 
terms being discarded. 

An algorithm can be derived from Eq. (37) which possesses 
an accurate estimate of the discretization error. Conversely, 
with a pre-defined level of acceptable discretization error 
this algorithm permits the selection of the load increment 
such that this level will not be exceeded. Thus, the in- 
cremental analysis proceeds with variable steps throughout 
the load-displacement history. 

The algorithm in question is based on a method of solu- 
tion of first-order differential equations attributed to 
Bulirsh and Stoer (Refs. 16,17). Two separate estimates of 
displacements are given at every other load increment. 

These values are averaged before the algorithm is repeated 
with a new starting point. The basic steps in the method 
are listed below and shown in Figure 9. 



30 


(1) U i+1 ) " {V ♦ dn [^up ]' 1 {P} 

(2) uA + 2 ) = up + 2d n [K T CA i+1 ) ] _1 (P) 

(38) 

(3) U? +2 } = (A ltl ) ♦ dn [K T CA iH . 2 3 ] - 1 {P) 

(4) U it2 > = \ (A^ 2 ) ♦ \ (A? t2 ) 

To estimate the discretization error o£ the above 
procedure, a simple cubic function is used as a test case. 
Consider the differential equation 

/ 2 

y = a + 2bx + 3cx 

which has the obvious solution 

, 2 3 

y = ax + bx + cx 

when y(0) = 0 

Applying steps (1) thru (4) when the step size is h, 

Y 1 58 ah 

y^ = 2ah + 4bh 2 + 6ch 3 
y 2 - ah + (ah + 4bh 2 + 12ch 3 ) 

Y 2 “ 2ah + 4bh 2 + 9ch 3 

6 = y 2 ~ y 2 = 6ch3 
But y(2h) exact = 2ah + 4bh 2 + 8ch 3 

3 3 

The error, e = ch * 0(h ) 


or expressed in terms of the difference, 6 between the two 
estimates of y 2 , 
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Extending this finding to the multidegree of freedom system, 


The error is defined as 

e - Max 


a A a B 
& i*2 , j ~ i+2 , j 

A i*2J 


(39) 


where j denotes the freedom number. 

The fact that the local discretization error e is 0(dn 3 ) com- 
bined with Eq. (39) can be used to either increase or decrease 
the load step size. Assuming that the predefined error toler- 
ance is denoted by TOL, then if 

e > TOL 

The step size is halved, and the procedure initiated from 
{A^} again. If 

8e < TOL 

the step size is doubled. 

A simple geomet rically nonlinear problem has been solved 
by the procedure indicated. The test case (Fig. 10) is a shal- 
low truss for which there is an analytical solution as well as 
other numerical solutions (Ref. 18). For comparison purposes 
the results from the modified incremental scheme are also dis- 
played in Figure 10. The number of times the stiffness matrix 
had to be inverted is approximately the same for both methods, 
and thus the cost of the solutions are comparable. The value 
of TOL was set at 0.1 for this example. The points A-D illus- 
trate how this tolerance was applied. When the load level 
reached point A the two subsequent increments produced the 
values indicated by points B and C. Since these failed to 
meet the error tolerance they were discarded, the increment 
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size was halved, and the satisfactory point D was obtained 
in the succeeding calculation. Although the present approach 
does not stay as close to the correct solution as the incre- 
mental method for most of the path, neither does it drift as 
far away at certain points. 
d. Overall Strategy 

In the present study the approach used for nonlinear 
analysis is the Newton-Raphson method. This method accommo- 
dates the modified incremental approach simply by prescribing 
a very large convergence tolerance factor. By this artifice 
only one iteration is performed at each load level. It is 
noted in Section 6b above that the modified incremental method 
is equivalent to one iteration of the Newton-Raphson method. 

Although no results for shell problems can be quoted for 
the initial value approach mentioned in the last section it 
is felt that this method holds great promise. This is particul- 
arly true for problems for which there are sudden changes in 
stiffness of the structure (e.g. an imperfect cylinder). It 
is hoped that this approach can be included in the same com- 
puter program as the Newton-Raphson method. This would facil- 
itate the choosing of an overall strategy by the user. 

7. CALCULATION OF LIMIT POINT 

As indicated previously, the determinant of the tangent 
stiffness matrix is zero at the limit point and consequently 
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the system of incremental equations does not possess an 
unique solution. This can be overcome by constraining the 
system in some way, such as through enforcement of a pre- 
scribed value of a particular displacement. This approach 
has been adopted in a variety of different forms by many in- 
vestigators in geometrically nonlinear finite element analy- 
sis, e.g., Refs. 19-23. The approach adopted here represents 
an extension of a scheme described by Zienkiewicz (Ref. 23). 

First, a representative degree-of -freedom (A^) is identi- 
fied. The central displacement of a plate or shell is an ap- 
propriate choice, the selection being input directly by the 
user. A change (6A^) in the magnitude of this displacement 
component is likewise specified. The value of the correspond- 
ing change in load intensity (<5n) is sought, this being equal 
to the difference in intensities between the (i)th and (i+l)th 
load intensities, i.e., 6p = - n^. The total incremental 

loading is given by 6n {P} . Three loading and support condi- 
tions are analyzed: (see Fig. 12) 

(1) No loading except for the force necessary to produce 

the prescribed value <SA . P * is calculated in the manner 

r r 

of a support reaction in this analysis. The calculated 

displacements are designated as {A'*'}. (This vector in- 
cludes the specified value SA r ) . 

(2) Application of the normalized load vector {F} with the 

ith degree-of-freedom held fixed (A^ = 0) . Calculate the 

(o') 

support reaction P r v 1 in this analysis and the displace- 
2 

ments (A }. 
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(3) The out of balance load vector (f), calculated on the 

basis of the displacements at the last load level (n^) 

is applied with the degree-of -freedom A^ held fixed. 

r 3 ^ 

Calculate the support reaction P^; J in this analysis 

and the displacements {A^}. 

In combining these cases it is useful to consider two 

separate conditions. The first combination considers cases 

(1) and (2) and excludes the effect of out of balance forces 

(Case 3) . For this circumstance the resulting displacement 

state is given by the displacements {A 1 } less the normalized 

displacements times the force P^ 1 ^ , i.e. 

P r p (l) 

{A 1 ) - -fo (A 2 ) 
r 


This gives the displacement state corresponding to zero reac- 
tive force at A^ for simple incrementation (no out of balance 
force correction) . 

The second combined condition takes into account the in- 
fluence of the out of balance force correction. By reasoning 
similar to that given above for the first combined condition, 
the resulting displacements are 

p( ' 3 ' 1 i * 

dn u > * <* > 


r 


For the modified incremental approach the out of balance 
loads have to be added over and above the loads of case (2) . 

It is important to keep these load vectors separate during the 
solution stage because at that point the reactions P^ 1 ^ , P^ 2 ^ 
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and P v 1 , and hence the factors of the above combinations, 
r 

are still unknown. 

Unfortunately, the combination of cases (1) and (3) 
gives no control over the displacement at the degree of free- 
dom A_ but in fact defines the value of the load there. In 
r 

consequence, at the limit point, where becomes zero 

(this is instantaneously the snapping through circumstance 
and the specified value 6A r is accomplished with theoretical- 
ly no required force p£^), the displacements become infinite. 
The problem is circumvented by neglecting the out of balance 
loads in the proximity of the limit point. The criterion 
used in deciding when to neglect the out of balance loads is 
the following 

P (3) > p(l) 
r — r 

Thus, in summary, the following combinations of cases are 

used. 

For P^ 3 ) < P^ 1 ) 
r r 


P (3) p(l) 

(1 + -X^) {A 1 } - {A 2 } + {A 3 } = { 6 A } (40) 

P r P r 


where {6A } is the increment of the total displacement state 

due to 6A - and for P^ 3 ^ > 

r ' r — r 


pd) 

(A 1 ) - -^y {A 2 } = { 6 A ) 


( 41 ) 


The change in load intensity can be determined after 


calculation of the increment of displacements as 
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6ti 


pU) 

r 


(42) 


or, since 


i + 1 



+ n i 


(43) 


It should be noted that after the limit point the ratio 
P^/P^ becomes negative, thereby reducing the load level 
until a minimum is reached at a snapback load, after which 
it increases again. 

To apply the prescribed displacements to the linear in- 
cremental equations the following approach (Ref. 24) is 
adopted. 

1) The origin of the displacement to be prescribed is 

changed. Thus, (6A r + A r ) replaces A r . By substitution 
in the equilibrium equations 


[ K T ] { A > = {P> - [K t ]{A} (44) 

where (A) consists of zeros everywhere except for freedom 

r which contains 6A . The modification to the load vector 

r 

{P} need only be carried out at freedoms coupled to the 
constrained freedom. 

2) A r is set to zero. This is accomplished by uncoupling 

the rth equation from the rest, by multiplying the diagonal 

20 

term by a large number, say 10 

Finally, after solution of the equations the reaction 
at freedom r is 

P r ' - lo2 ° A r K T 


(45) 
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This procedure is obviously valid if there is more than 
one prescribed displacement. 

8. CALCULATION OF BIFURCATION STATE 
a. Linear 

The method of calculating the intensity of load (n c ) for 
bifurcation and the associated mode shape is based on the 
condition of zero value of the determinant of the tangent 
stiffness at this point. This condition holds because of 
the lack of uniqueness of the load displacement path at the 
bifurcation point, resulting in the singularity of [K,p] . The 
rank of the singularity depends on the number of paths present 
at the bifurcation point. In this development attention is 
restricted to cases where only two such paths are present. 

For these cases the singularity is of rank one. 

It should be noted that past reports of this research 
(Refs. 9 and 10) have emphasized approaches to calculation of 
bifurcation which depend on the condition stated above. The 
following description differs from these principally in the 
computational mechanics. 

The problem of linear stability, characterized by a 
linear pre-buckling state, is examined first. The condition 
cited above can be written for the general nonlinear represen- 
tation as 

|K t | = | [K o ] + [Kj] + [K 2 ] | = 0 (46) 

where | | denotes the determinant of the indicated matrix. 

To transform this to the linear buckling problem the matrix 
[K 2 3 is neglected. This matrix is in fact zero in the absence 
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of coupling between inplane and out -of -plane displacements 
on the fundamental path, as in the case of flat plates. 
Secondly, the [K^] geometric stiffness matrix is assumed to 
vary linearly with the load level n. This assumption is 
justified if the distribution of internal forces in the 
structure does not change along the fundamental path. 

The problem therefore becomes one of finding the critical 
load level (n ) such that 


|[K 0 ] + n c [^(4)] I = 0 


(47) 


where {A} are the displacements corresponding to a normalized 
load vector (n = 1) . 

The determinant expressed in Eq. (47) is first established 
at two points, n = 0 and n = 1. These two points furnish suf- 
ficient information for starting a Lagrangian interpolat ive 
scheme which subsequently progresses in an iterative fashion. 

An analytical expression of polynomial form is sought for n 
in terms of already-calculated values of load intensity (r^) 
and the determinants (Det.j) expressed by Eq. (47) at these 
intensities. This form is obtained by Lagrangian interpola- 
tion and writing the result as an equation which gives the 
intensity (n^ +1 ) for a zero value of the corresponding 
determinant one has 


i + 1 


l 

I 

j=l 


i(Mj) 

( n 

k=l 


Pet, k 
Det . k-Det . j 


■) 


(48) 


The process is shown schematically in Fig. 12. Each iteration 



39 


consists of the following steps 

1) Scale element geometric stiffnesses [K^] by and 

add into total stiffness matrix which is initialized 
to [K q ]. 

2) Decompose total stiffness matrix into upper triangular 
form by Gaussian Elimination procedure. 

3) Evaluate determinant by multiplying the diagonal terms. 

4) Perform Lagrangian interpolation to obtain a better 
estimate of the load required to yield zero determinant. 

The iteration procedure continues until two successive values 
of n are within a predefined tolerance. 

After convergence of the above procedure, the eigen- 
vector can be calculated with very little extra computational 
effort. The problem formulation is now that of finding a 
displacement vector { A > that satisfies the following relation- 
ship . 

[K o + n c K i CAD ] { A> = 0 (49) 

The decomposition of these equations into upper triangular 
form has already been performed during the last iteration of 
the eigenvalue procedure. To insure a nontrivial solution 
for {A), before proceeding with the backsubst itution , the 
right hand side is initialized to zero except for the last 
entry which is given a value of unity, as is the last diagonal 
term of the decomposed matrix. By this artifice the last 
displacement in the structure is forced to take a value of 
unity, and thus the eigenvector is normalized with respect 
to this displacement. Note that this method will yield a 
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nontrivial eigenvector flat plates and beams, only if the 
last freedom in the structure is a bending displacement. 

The user must ensure this by the way he numbers the nodal 
points in the idealization. 

The choice of the last equation for discarding is 
rather arbitrary. In the present case the choice is dic- 
tated by computational simplicity. Wilkinson (Ref. 25) has 
shown that for some matrices inaccurate results are obtained 
unless one particular equation is discarded. Because of 
this he suggests a symmetric strategy: inverse iteration. 

In this method an approximate eigenvector {A^} yields an 
improved approximation {A^^} by the following procedure. 

[K o + n c K l (A)] {A i+l } = {A i } (50:) 

As already noted, the matrix on the left hand side of Eq. 

(SO) has been triangularized during the last iteration of 
the eigenvalue procedure. Thus, successive resolutions of 
the equations indicated by Eq. (50) are called for with the 
right hand side consisting of the displacement vector from 
the last iteration. This procedure is computationally quite 
inexpensive because the elimination procedure has to be per- 
formed only on the right hand side. This algorithm has been 
found (Ref. 26) to converge in a few cycles in most cases. 

In the present study the eigenvector obtained by dis- 
carding the last equation has proven to be highly accurate. 

If the accuracy of this solution proves to be suspect for a 
class of problems not yet solved the inverse iteration scheme 
described above can be incorporated with a minimal effort. 
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b. Nonlinear 

The approach of the last section, wherein the determin- 
ant of the tangent stiffness matrix [K T ] is equated to zero, 
is again used as the basis for predicting nonlinear buckling. 
In the presence of nonlinearities in behavior along the 
fundamental path the simplifying assumptions exploited to 
simplify the problem to a linear one are no longer valid. 

In other words, [K 2 ] is retained, and [K-^j is not assumed 
to be linearly dependent on n. Consequently, the nonlinear 
buckling problem is far more complex to solve than its 
linear counterpart. 

In the present study, the load-deflection curve of the 
structure is followed up to bifurcation. An incremental 
procedure is used for this, so that the determinant of the 
tangent stiffness matrix can be monitored during the analy- 
sis. Constant load increments are used throughout, so that 
Lagrangian interpolation with equal intervals is used to 
predict where the determinant goes to zero. A comparison 
of linear and nonlinear buckling with respect to the deter- 
minant-load level plot, is shown in Figure 12. When the 
load level has reached a point such that buckling will occur 
during the next load increment a linear buckling analysis is 
performed. The last tangent stiffness matrix calculated 
takes the place of the linear stiffness matrix in Equation 
( 49 ). 
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9. NUMERICAL ANALYSES 

Numerical solutions to four nonlinear problems are 
described in this section. These have been chosen because 
of the availability of comparison solutions. 

The first problem is the case of a constant thickness 
(t) isotropic square plate subjected to uniform uniaxial 
compression (see inset. Fig. 13). The boundary conditions 
of the finite element analyses consist of simple support 
with respect to flexure (i.e., zero normal displacement 
and normal bending moment along the edges) and complete 
freedom for the edges to displace in the plane of the plate. 
Initial imperfections are present as represented by a nor- 
mal displacement of the middle surface, distributed in 
double sinusoidal form, with a central amplitude O.lt. Be- 
cause of the symmetry of load, geometry and structural response 
the analyses are performed for only a quadrant of the plate. 

The solutions of the problem are presented in Fig. 13 
in the form of a plot of the ratio of total applied load (P) 
and the critical load for linear stability (P ) versus the 
ratio of the central deflection (w) and the plate thickness. 

It is seen that the range of solutions extends to P/P cr = 2. 

A classical solution for this problem due to Coan (Ref. 

27) is presented along with the finite element solutions, that 
are in each case based upon load increments of P cr /10. It 
is appropriate to first consider the results obtained by use 
of a full Newton-Raphson iterative scheme, in which itera- 
tions in each load increment are continued until convergence 
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is achieved. These are given for idealizations consisting 
of two elements (the lxl grid in the quadrant, consisting 
of 2 triangles) and 8 elements (the 2x2 grid) . It is seen 
that these converge to a solution to the right of the clas- 
sical solution, i.e., a solution that is less stiff than the 
classical solution. The modified incremental approach (the 
tangent stiffness procedure combined with a one-step Newton- 
Raphson correction) is applied to the 8 -element representation. 
It is interesting to note that for the chosen load increment 
the difference between the full Newton -Raphson and the less 
expensive modified incremental results is negligible. 

The discrepancy between the classical solution and the 
above finite element solutions is due in part to the constraint 
of the former to linear inplane displacement of the loaded 
edges. This gives a stiffening effect. It is not customary 
in finite element analysis to enforce such an idealized state 
upon the edge conditions and in the above-described finite 
element solutions the inplane edge displacements were not con- 
strained to give a linear variation. When this constraint is 
enforced for a 2-element (lxl mesh) finite element solution, 
stiffer results than Coan’s are obtained (see values denoted 
by x in Fig. 13), as expected. 

The second nonlinear analysis problem is shown in Figure 
14. A cylindrical segment is fixed along all sides and is 
subjected to pressure loading. The solutions to be discussed 
involve the calculation of the central deflection (w c ) as a 
function of the applied pressure. Brebbia and Connor (Ref. 28) 
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first solved this problem with use of rectangular shell elements 
based on assumed displacements, i.e., a stiffness formulation. 

Lien (Ref. 29) also employed rectangular shell elements. Prato 
(Ref. 30) applied a mixed, Reissner-type shell element to this 
problem. Dhatt (Ref. 31) utilized a discrete -Kirchoff approach. 
Finally, it should be noted that Connor and Morin (Ref. 32) have 
reported numerical results for this problem. 

Results from the present formulation, for a 3 x 3 (18-element) 
gridwork, are shown in Fig. 14. Advantage is taken of symmetry, 
with a 3 x 3 (18 -element) gridwork in a quadrant of the shell. 

None of the alternative numerical solutions from Refs. 28-32 
gives tabulated data and it is felt that the graphical represen- 
tations in these references, except for Ref. 31, cannot be 
transcribed with the accuracy associated with plotted comparisons. 

The basic difference in these results, up to a pressure of 
approximately 0.150 is accounted for by the difference in linear 
solutions, the Ref. 31 results being based on a different varia- 
tional principle with 9-term expansions for u, v, and w and 
independent quadratic expansions for the middle surface slopes. 
These same considerations may have influence on the nonlinear 
solutions. Although, as noted above. Ref. 32 does not give data 
that is accurate enough for plotted comparisons , it appears to 
the authors to have a solution that corresponds more closely to 
the present work. 

The third nonlinear analysis problem pertains to the 
spherical cap shown in Fig. 15, which is subjected to a concen- 
trated load P at the crown. The support conditions involve the 
restraining of displacements normal to the edge and normal to 
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the shell itself. Thus translations along the edge of the 
shell is permitted. This type of support condition is some- 
times described as ’hinged’. 

This structure exhibits a snap -through behavior under the 
indicated load with eventual recovery of stiffness at higher 
displacements. A classical (series) solution for the problem 
has been obtained by Leicester (Ref. 33) and an alternative 
finite element solution by Dhatt (Ref. 31) who, as was noted 
previously adopts a discrete -Kirchoff approach. 

Again, a 3 x 3 grid (18 -elements) is employed within a 
quadrant of the structure. Results are shown in Fig. 15 for 
the central deflection as a function of the applied load. The 
agreement among the three solutions is remarkably good through 
two limit points, up to where a stiffening behavior is again 
experienced. As is usual with the incremental approach the 
displacements tend to lag behind the true solution for increas- 
ing loads. Despite this, the value of the snap-through load 
is predicted very accurately by the present method. 

To verify the accuracy of the quadrilateral element 
formulation the problem of a square clamped plate under pressure 
loading is considered. A 4 x 4 idealization is used for one 
quarter of the plate as shown in Fig. 16. The analysis proceeds 
to a load level that produces considerable nonlinear behavior. 
The modified incremental approach is used to solve the nonlinear 
equations. The results are displayed in Fig. 16 along with the 
classical solution of Levy (Ref. 35) , and other finite element 
solutions. Again the present results lag slightly behind the 
true solution. 
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10. CONCLUDING REMARKS 

This report has presented the formulation, for geometrical- 
ly nonlinear analysis, of a curved triangular thin shell finite 
element together with computational approaches which enable 
the use of the element formulation in a variety of thin shell 
instability calculations. The linear aspects of the element 
formulation are described in a separate report (Ref. 1). The 
types of instability calculations which were effected included 
bifurcation and limit point analyses. 

The basis of the finite element formulation was a cubic 
polynomial field within the element. Although this is close 
to the simplest formulation that can be attempted for triangular 
elements which evidence bending, it is not efficient when used 
without modification due to the discontinuity of normal slope 
across element interfaces. The requirement of slope continuity 
is therefore written as a constraint condition and is appended 
to the global equations by means of the Lagrange multiplier 
technique. Numerical analyses presented in this report and in 
Ref. 1 show that efficient solutions are obtained in this manner 
for an extensive range of comparison problems. Displacement 
fields which were simpler than those used in the linear formula- 
tion were used in the development of terms for geometrically 
nonlinear behavior. This approach, known as an ''inconsistent" 
formulation, did not significantly prejudice the accuracy of 
the problems solved numerically. 

The cornerstone of all of the types of nonlinear analyses 
presented herein was the modified incremental method, which 
combines a conventional incremental, or tangent, step with a 
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one-cycle Newton-Raphson calculation. Numerical comparisons 
showed that the adopted approach strikes a balance between ac- 
curacy and computational efficiency for the problems solved. 

An alternative approach was presented for the solution of the 
nonlinear algebraic equations, based upon treatment of these 
equations as an initial value problem. This approach was 
demonstrated for a simple test problem but was not applied to 
the large-scale numerical analyses. It features the rational 
selection of step size and the attainment of solution accuracy 
within prescribed limits and therefore appears worthy of further 
study and improvement . 
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Figure 1. Forms of Nonlinear Force -Displacement Response 
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Fundamental Path 

Imperfect Structure 

Response 



a. Response of Actual b. Limit Point Phenomenon 

Structure Insensitive 
to Imperfections 


Figure 2. Critical Load Phenomena 



c. Cylinder in compression 


d. Shallow truss 


Figure 3 


Examples of Structures Sustaining Various Forms of 
Instability and Nonlinear Force-Displacement 









Figure 4. Triangular Thin Shell Element 


Figure 5. Axial Member-Geometrically Nonlinear Behavior 



O Retained nodal point 
X Eliminated nodal point 

Figure 6. Quadrilateral for Inconsistent Formulation 
of Nonlinear Matrices 








Deflection at vertical load (inches) 


Figure 10. Solution Comparisons -Shallow Truss Problem 





59 



'Out of Balance' Loads 


Figure 11. Method of Limit Point Calculation 



Determinant of 


60 




61 



Figure 13. Imperfect Plate Under Edge Compression 
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Figure 15. Spherical Cap Under Concentrated Load at Crown 
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Figure 16. Clamped Square Plate Under Pressure Load 

NASA- Langley, 1975 CR-2483 



